home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
toolbox
/
tfd
< prev
next >
Wrap
Text File
|
1996-04-09
|
7KB
|
251 lines
# function [y,min_t,max_t]=tfd(x,fa,sigma,threshold)
#TFD Choi-Williams- and Wigner-Ville time-frequency distributions.
#
# USAGE: y=tfd(x,fa,sigma,threshold)
# y=tfd(x,fa,sigma)
# y=tfd(x,fa)
# y=tfd(x)
# [y,min_t,max_t]=tfd(...)
#
# The default values for the non-defined parameters are
#
# fa = round(length(x)/2)
# sigma = inf (i.e. a Wigner-Ville distribution)
# threshold = 0.001
#
#
# INPUT:
# x Complex input signal (column- or row vector). If signal is
# real, use hilbert(x).
# fa Length of frequency axis. (Default value: half the
# signal length.)
# sigma Degree of smoothing. The _smaller_ sigma, the _more_
# the crossterms are suppressed in the Choi-Williams
# distribution. Nice value to start with: sigma=1.
# If sigma increases the result tends to a Wigner-Ville
# distribution. If sigma=inf (default) then an exact
# Wigner-Ville distribution is returned.
# threshold Minimum weightfactor (>0). In the Choi-Williams tfd the
# exponential weighting function in expands to
# infinity. However terms that are multiplied with a very low
# weighting do add to the computation time, but do not add
# significantly to the final result. Using this parameter,
# only terms that have a weighting larger than "threshold" add
# to the result.
#
# OUTPUT:
# y - Contains the distribution. Each row represents
# a frequency, each column a time instant.
# min_t - first time-instant of distribution
# max_t - final time-instant of distribution
#
#
# EXAMPLE 1:
# x=exp(2*pi*j* 6*((1:40)/40).^2); # Make a chirp
# fa=20; sigma=inf;
# [y min_t max_t]=tfd(x,fa,sigma); # Make a WVD
# contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of WVD
# hold on
# plot((min_t:max_t),12*(min_t:max_t)/40) # Plot instantaneous frequency
# hold off
#
# EXAMPLE 2:
# figure(1)
# x=exp(2*pi*j* 6*((1:80)/80).^2)+... # Two 'crossing' chirps
# exp(2*pi*j*(1:80)/80.*(50-30*(1:80)/80));
# fa=40; sigma=inf;
# [y min_t max_t]=tfd(x,fa,sigma); # Make a WVD
# contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of WVD
# figure(2);sigma=1;
# [y min_t max_t]=tfd(x,fa,sigma); # Make a CWD
# contour((min_t:max_t),(1:fa), flipud(y))# Contour plot of CWD
# hold on
# plot((min_t:max_t),12*(min_t:max_t)/80) # Plot instantaneous frequencies
# plot((min_t:max_t),(50-60*(min_t:max_t)/80))
# hold off
#
# Author: Rene van der Heiden, R.vanderHeiden@fel.tno.nl
# Ref: - H.I. Choi and W.J. Williams; IEEE Trans. Acoust., Speech, Sign.
# Proc., Vol. ASSP-37, pp 862-871, June 1989
#
# Translated to Rlab, 11/11/95, Ian Searle
#
# This software may be freely used and modified for research and development
# purposes. If you wish to use it for commercial gain please contact me.
# I provide ABSOLUTELY NO WARRANTY for this software.
#
# Copyright: Rene van der Heiden
#
tfd = function ( x, fa, sigma, threshold )
{
local (x, fa, sigma, threshold)
#
# Transpose signal if it is a column:
#
NROW = x.nr; NCOL = x.nc;
if (NCOL == 1)
{
x=x.';
LENGTH=NROW;
else
LENGTH=NCOL;
}
#
# Input check and, if necessary, assignment of default values
#
def_sigma = inf();
def_fa = floor(LENGTH/2);
def_threshold = 0.001;
if (nargs == 1)
{
sigma = def_sigma;
fa = def_fa;
threshold = def_threshold;
else
if (nargs == 2)
{
sigma = def_sigma;
threshold = def_threshold;
else
if (nargs == 3)
{
threshold = def_threshold;
}
}
}
if (threshold <= 0)
{
disp("Parameter \"threshold\" is chosen less than or equal to zero. ");
disp("Using threshold="+num2str(def_threshold)+" instead.");
threshold = def_threshold;
}
# In the next line we determine the size of the
# frequency window we will really use: window_use. (it should be uneven)
#
# eg. if fa equals 8 then window_use becomes 7 and
# if fa equals 7 then window_use also becomes 7.
#
window_use = 2*floor((fa-1)/2)+1;
#
# Calculate the number of timepoints at which the distribution will be
# evaluated:
#
length_cwd = LENGTH - window_use + 1;
#
# Initialize cwd:
#
cwd = zeros(length_cwd,fa);
#
# Calculate the maximum size of the windows over which we smooth:
#
wm = 2*floor((window_use-1)*sqrt(-log(threshold)/sigma));
#
# Calculate for each tau and mu what the weights are:
#
wgts = zeros((window_use-1)/2+1,wm+1);
mu = -wm/2:wm/2;
for (tau in 0:(window_use-1)/2)
{
if ((tau != 0) && (sigma != inf()))
{
wgts[tau+1;1:wm+1] = exp(-mu.*mu*sigma/(4*tau*tau));
else
wgts[tau+1;1:wm+1] = (mu < 1);
}
}
#
# n is the time variable:
#
min_t = (window_use+1)/2;
max_t = LENGTH - (window_use-1)/2;
for (n in min_t:max_t)
{
#
# tau runs over the right part of window_use:
#
for (tau in 0:(window_use-1)/2)
{
#
# Summation over the second window to smooth the distribution:
#
wm_tau = 2*floor(2*tau*sqrt(-log(threshold)/sigma));
mu_begin = -min(wm_tau/2,-1+ min(n+tau,n-tau));
mu_end = min(wm_tau/2,LENGTH-max(n+tau,n-tau));
weights = wgts[tau+1;mu_begin+wm/2+1:mu_end+wm/2+1];
#
# Scale the weights:
#
weights = weights/sum(weights);
#
# And calculate the distribution:
#
x1 = x[n+tau+mu_begin:n+tau+mu_end];
x2 = x[n-tau+mu_begin:n-tau+mu_end];
cwd[n-(window_use+1)/2+1;tau+1] = sum(weights.*x1.*conj(x2));
}
}
#
# Because the numbers in the CWD must be hermitic as a function
# of TAU, we can conjugate the found values of the CWD and copy
# them to the other half of the CWD.
#
if (rem(fa,2) == 0)
{
#
# if fa is even then write a zero in the cwd:
#
cwd[n-(window_use+1)/2+1;fa/2+1] = 0;
for (tau in fa/2+2:fa)
{
cwd[;tau] = conj(cwd[;fa+2-tau]);
}
else
for (tau in (fa+1)/2+1:fa)
{
cwd[;tau] = conj(cwd[;fa+2-tau]);
}
}
#
# And calculate the one dimensional simultaneous FT:
#
y = real(fft(cwd',fa));
return << y=y; min_t=min_t; max_t=max_t >>;
};